pkgs <-c('tidyverse','lintr','sf','raster','viridis','cowplot','ggpubr','ggspatial','haven')
lapply(pkgs, function(x){library(x, character.only = T)}) 

setwd('D:\\WorkingPaper-Series\\Township_Light\\Replication_Files\\workdata')

rm(list = ls())

town_light_tbv <- read_dta('township_panel_final.dta')

county_tbv <- town_light_tbv %>% 
  group_by(idcode) %>% 
  mutate(n_tbv = max(N_village)) %>% 
  distinct(idcode, .keep_all = T) %>%
  group_by(countycode) %>%
  summarise(n_tbv = sum(n_tbv))

# load in the maps
county_shp <- read_sf("D:\\WorkingPaper-Series\\Township_Light\\Replication_Files\\workdata\\county.shp")
county_shp <- county_shp %>% dplyr::select(PAC) %>% rename(countycode = PAC) %>% 
  left_join(county_tbv, by = 'countycode')

china_shp <- read_sf("D:\\WorkingPaper-Series\\Township_Light\\Replication_Files\\workdata\\nation_line.shp")
china_shp <- china_shp %>% dplyr::select(geometry) %>%
  st_transform(crs = st_crs(county_shp))


# plot
breaks <- c(1,5,10,20,50,100)

main_map <- county_shp %>% 
  mutate(fire_discrete = as.factor(findInterval(n_tbv, vec = breaks))) %>%
  ggplot() +
  geom_sf(aes(fill = fire_discrete), color = "grey50", size = 0.1) + 
  coord_sf(ylim = c(15, 55)) +
  theme_bw() +
  scale_fill_manual(
    values = c("0" = "gray95", 
               "1" = "#F7A300",
               "2" = "#F77B00", 
               "3" = "#F74B00", 
               "4" = "#C41C00",
               "5" = "#7F0000", 
               "6" = "#4D0000"), 
    labels = c("0" = "Counties with 0 TBV",  
               "1" = "1 - 5 TBVs",
               "2" = "5 - 10 TBVs",
               "3" = "10 - 20 TBVs",
               "4" = "20 - 50 TBVs",
               "5" = "50 - 100 TBVs",
               "6" = "TBVs Above 100"),
    na.value = "white",
    name = NULL) +  
  theme(
    legend.position = c(0.1,0.15),
    legend.background = element_rect(color = 'black', size = 0.1,linewidth = 0.1, fill = 'white'),
    text = element_text(family = 'serif', size = 10),
    legend.text = element_text(family = 'serif', size = 10),
    axis.text =  element_text(family = 'serif', size = 10),
    panel.grid.major = element_line(color = gray(0.5), 
                                    linetype = 'solid', 
                                    linewidth = 0.1)
  )


south_sea_map <- china_shp %>% 
  ggplot() +
  geom_sf(fill = "gray95", color = "grey70", size = 0.05) + 
  coord_sf(
    xlim = c(105, 125),
    ylim = c(2, 25),
    expand = FALSE
  ) +
  theme_void() +
  theme(
    panel.background = element_rect(fill = "white", color = "black", size = 0.5),
    panel.border = element_rect(color = "black", fill = NA, size = 0.5)
  )

final_map <- ggdraw() +
  draw_plot(main_map) +
  draw_plot(
    south_sea_map,
    x = 0.8,      
    y = 0.05,      
    width = 0.2,   
    height = 0.25   
  )


